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Abstract: When using a formulation of Smooth Particle 
Hydrodynamics (SPH) which conserves momentum exactly the 
motion of the particles is observed to be unstable to negative 
stress. It is also found that under normal circumstances a 
lattice of SPH particles is potentially unstable to transverse 
waves. This document is a summary of a detailed report (Morris 
1994) investigating the nature of these and other instabilities 
in depth. Approaches which may be used to eliminate these 
instabilities are suggested. It is found that the stability properties 
of SPH in general improve as higher order spline interpolants, 
approximating a Gaussian, are used as kernels. 

1. Introduction 

The early applications of SPH, Lucy 1977 and Gingold and 
Monaghan 1977, were to problems involving a compressible gas 
which always had a positive gas pressure. The equations govern- 
ing the motion of the particles, when written in a form which 
conserves momentum exactly, result in particles repelling each 
other with equal and opposite forces. As particles approach, the 
density increases, the pressure increases, and the particles tend 
to repel each other. When applied to different problems where 
the stress can become negative the momentum conserving form 
of SPH is observed to become unstable to short wavelength per- 
turbations. For negative stress, the particles no longer repel, 
but attract. Each particle is, in effect, at the bottom of a po- 
tential well, and it becomes possible for the particles to pair 
up and "slide" into each others wells, causing "clumping" ini- 
tially and subsequently disrupting the solution. The nature of 
this instability is quite different from instabilities observed in ex- 
plicit finite difference techniques when too large a time step is 
used. This SPH instability is a consequence of the particles be- 
ing free to move under negative stress and is present even if the 
time integration is exact. Additional stability conditions, such 
as the Courant Friedrichs Levy condition must also be respected. 
Recently people have commented upon this instability in con- 
nection with many problems (e.g. elastic now in Swegle 1992), 
but much earlier Phillips and Monaghan 1985 had detected and 
studied the problem in connection with magnetohydrodynamics 
(MHD). Very little analytical work has been done in the study of 
the stability properties of SPH and previous methods (see Mon- 
aghan 1989) actually break down in the region where the stability 
we seek to understand first appears. The report (Morris 1994) 
extends the stability analysis of SPH in order to not only under- 
stand the nature of potential instabilities but to give insight into 
how methods giving most accurate results may be formulated. 

2. Seeking the Source of the Instability 

In order to understand the problem better we must seek the 
fundamental source of instability. We have seen that the reported 
instability occurs when a stress tensor is implemented in a form 
which conserves momentum exactly. We will be considering one 
dimensional flows initially, with symmetric kernels. Typically, we 
might use a Gaussian kernel, 

or the cubic spline interpolant, 
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Choosing a symmetric kernel guarantees exact conservation 
of momentum, since inter-particle forces then form action- reaction 
pairs. 

Let us consider the specific case of one-dimensional MHD 
flow, with a constant magnetic induction Bq parallel to the 
direction of the flow, say along the a>axis. Since the magnetic 
induction is constant, the equations of MHD reduce to those for 
hydrodynamics. However, the corresponding SPH equations may 
be written, 

dv a ^ I 'Pa - 2^B<> Pb - ^;B<>\ dWab 

Pa = ^2rn b W ab , p a = c 2 p a , 
b 

W ab = W(x a -x b ,h). 

So, in this case, the SPH equations for magnetohydrodynamical 
equations are equivalent to those for pure hydrodynamical flow 
with the pressure adjusted by a constant. Similar analysis may 
be done for any problem involving a stress tensor of similar 
form. So, the stability properties of one-dimensional Smooth 
Particle Magnetohydrodynamics (SPMHD) will be exhibited by 
the equations, 




Pa = ^2 r m b W ab , p a = C 2 p a + P, 



b 

where P = —Bq/2/io. Similar equations, with a different expres- 
sion for P will be obtained by considering the one-dimensional 
case of other applications of SPH. 

3. Stability Analysis and SPH 

A simple way to analyse the stability of an implementation 
of SPH is to consider the dispersion relation for linear, sound 
waves propagating in a one-dimensional flow. We can imagine 
an infinite line of identical particles oscillating about a constant 
mean separation. Such a wave may be described by, 

x a = a Ax + X exp (ikaAx — iwt) , (5) 

We now linearise Eqs. (4) by substituting Eq. (5) and neglecting 
all but first order terms to obtain, 




where we have chosen, P = c 2 po (Jft — 1). I have assumed W to 
be even (a more general expression is given in Morris 1994). This 
choice of P gives, 

Va = c 2 Po (jjL - 1 + 9t) . (7) 

We can think of the quantity c 2 podt as being the "background" 
pressure upon which the linear wave perturbations occur. Thus, 
if 9? > 0, then p a should be positive for small perturbations and 
if 9i < 0, then p a may become negative. !ft will take different 
values depending on the equation of state being considered. For 
standard SPH, $t = 1. For SPMHD, $t = 1 - B 2 /2n c 2 p . For 



incompressible SPH (see Monaghan 1994), ^ : 
summation formula, Eq. (6) may be written, 
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4. Stability of Another Formulation 

It is possible to develop an alternative formnlation to Eq. (4) 
which has stability properties independent of the background 
pressure, 

dv a ^ Pb ~ Pa dW ab 

= — > m.i. . 

dt ^ 



PaPb 



dx„. 



(9) 



Since the differences of the pressures are taken in the momentum 
equation, Eq. (9), the value of !ft is irrelevant. However, the 
forces exerted upon the particles according to Eq. (9) no longer 
form action-reaction pairs but are estimates of the local pressure 
gradient. Thus, this formulation does not conserve momentum 
exactly, except in the continuum limit (ie: h —> and infinitely 
many particles). This pressure difference formulation has a very 
simple equation governing its stability, 



dW 

smkAxj——(Axj,h) 



(10) 



5. The Numerical Sound Speed 

Ideally, we have w 2 = c 2 k 2 , thus it is useful to define, 



k 2 c 



where c is the analytic sound speed. Ideally, of course, C num 
should be as close to one as possible, so that sound waves 
are propagated realistically. If C 2 um becomes negative, the 
perturbations in the numerical solution are no longer travelling 
waves, but are exponentially growing and decaying disturbances. 

6. Results 

It is clear that the dispersion relation is periodic in k with 
period 27r j Ax. This is the result of aliasing, whereby waves of 
wavenumber k and k + 2it j Ax look identical at a series of points 
separated by Ax. It has been established experimentally that the 
"problem" wavenumber is 7r / Ax (a wavelength of two particle 
spacings). For k = 7r/Ax, Eq. (8), becomes, 
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Eq. (11) shows that, for wavenumber tt/Ax, the sign of w 2 , and 
hence C 2 um , changes sign with Thus, for any symmetric 

kernel, as the background pressure c 2 podt changes sign, it will 
change stability. There can be no kernel which is stable for 
both signs of using this implementation of SPH. Thus, to 
achieve general stability for different background pressures, we 
must modify the kernel as !ft changes sign. Let us consider 



Eq. (11) for the simple case where W(Axj, h) = for \j\ > 2. 
Then, 

9 8mc 2 $d 2 W , A 
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So, for this case, the sign of d 2 W/dx 2 ( Ax, h) must change 
when the background pressure changes sign. This is particularly 
important to remember if we are seeking a kernel which changes 
continuously with to provide stability. The appropriate kernels 
for being slightly positive and slightly negative have different 
signs for d 2 W/dx 2 , so d 2 W/dx 2 must be zero at the nearest 
neighbour when = for continuity. In the case of the Gaussian, 
we find that w 2 = when h = 2~Ax. Thus, varying the 
smoothing length may be sufficient to ensure stability for this 
wavenumber. A useful solution, however, must be stable for all 
wavenumbers. 

We can see from the form of Eq. (6) that for any given h and 
k, C 2 um varies linearly with In general we wish C 2 um to vary 
only slowly with since, for the exact system, it is independent 
of It is informative to consider the quantity dC 2 um / ddt, since 
this reflects the degree of dependence of the stability properties 
on For sufficiently large !ft we have that C 2 um oc !ft where the 
constant of proportionality (for a given k and h) is dC 2 um / ddt. 

The limit of Eq. (8) as k approaches zero is of particular 
significance. It is in this limit that we expect our numerical 
method will resolve the wave best, since it corresponds to an 
infinite number of particles fitting into each wavelength. It can 
be shown that, in this limit, C 2 um is perturbed from 1 by a term 
which is linear in !ft and proportional to the Fourier transform 
of the kernel evaluated at 27r / Ax. The Fourier transform of the 
Gaussian kernel is, of course a Gaussian in kh, and, accordingly, 
falls off quite rapidly with kh. The Fourier transforms of kernels 
with compact support, however, do not fall off so rapidly, and we 
expect the sound speed to depend more strongly upon !ft for such 
kernels. Numerical results obtained using direct summation in 
Morris 1994 of Eq. (6) confirm that C num depends more strongly 
on !ft for spline interpolated kernels. In fact, there are choices 
of k, Ax and h for which dC 2 um /ddt is negative. Thus, a 
sufficiently large choice of !ft for such parameters results in C 2 um 
becoming negative, producing an instability. This result was not 
anticipated since it had been thought that SPH would be stable to 
any positive stresses. This instability resulting from over-pressure 
is of long wavelength and, correspondingly, can have a serious 
affect upon the accuracy of large scale phenomena. The use of 
smoother spline approximations to a Gaussian kernel increases 
the over-pressure required to induce the instability. More details 
on these results may be found in Morris 1994. 

7. Mutating the Kernel 

We have seen that no single kernel can provide us with the 
stability we seek. We need some way of modifying the present 
Gaussian shaped kernel. Any kernel, W(x, h) we choose must 
have the following properties (in one-dimension): 

/ + CO 
W(x, h)dx = 1, i!T W(x, h) = S(x), (13) 
-CO 

and, if we wish momentum to be conserved exactly, 



W(x,h) = W(-x,h). 
Suppose that our "favourite" kernel is, 

W 0{ x,K)= l f ^), 



(14) 



(15) 



We can now define an infinite series of symmetric kernels satisfy- 
ing Eqs. (13) of the following form, 



W n ( X ,h): 
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(16) 



where A n is chosen to satisfy Eq. (13). Now we have a set of 
kernels which we can combine to form a kernel with the necessary 
properties for stability, 



where, 



W(x,h) = J2B n W n (x,h), 



(17) 



(18) 



We need to develop some method for choosing the co- 
efficients B n snch that the resulting kernel has the desired 
stability properties. One possible approach is to choose several 
points that we want w(k) to pass through. We then solve for 
the same number of free co-efficients as we have conditions by 
linearising about an initial guess and solving the resultant linear 
system of equations. We continue this iterative process until 
sufficient convergence is achieved. There are some conditions 
which must be respected in order that the system converges. 
These, and other details, are explored in Morris 1994. 

8. The Pressure Differencing Formulation of SPH 

The form of Eq. (10) ensures that the formulation Eq. (9) 
is always stable. It can also be shown that the sound speed 
for this implementation depends less strongly upon the parti- 
cle spacing and smoothing length than the exactly momentum 
conserving formulation. This implementation, of course, has sta- 
bility independent of !ft since all pressures appear as differences in 
Eq. (9). However, it does not conserve momentum exactly and, 
thus, is not as good as the momentum conserving formulation 
(see Steinmetz and Mueller 1993) for modelling strong shocks. 
The pressure differencing formulation exhibits much larger os- 
cillations when modelling strong shocks, though these might be 
tamed by using a different form of viscosity. It is possible to 
devise a shock switch to detect and track shocks such that this 
viscosity is only present at the shock front. The exactly momen- 
tum conserving method calculates the net particle forces as the 
sum of individual repelling interactions. This is less subject to 
oscillations than a method where the forces exerted on the par- 
ticles are calculated from an estimate of the pressure gradient. 
When modelling stronger shocks, the pressure differencing formu- 
lation becomes less accurate as exact conservation of momentum 
becomes more necessary. 

9. Stability Analysis with Variable Particle Spacing 

In practice, a simulation will not involve particles oscillating 
about a constant separation. In order to deal with a number of 
different particle spacings we could solve for a kernel appropriate 
to each particle spacing and then use symmetric combinations of 
these when particles interact. If the expression Eq. (6) is used 
to find appropriate kernels, it is found that the numerical sound 
speed is significantly reduced. Effectively, the kernel is being 
changed continuously in such a fashion as to slightly oppose the 
propagation of the wave. If, however, we reformulate the stability 
analysis to include the assumption that the kernel is a function 
of the particle spacing, the sound speed will be more accurate. 
We need an appropriate variable to reflect the local spacing of 
particles. The distance to a nearest neighbour and other similar 
quantities are quite clumsy to use for such analysis. The specific 
volume, 

K = -, (19) 

Pa 

is much more suitable. We now consider Eqs. (4) with the 
consideration that, 



W ai = W(x ai , h, V ai ), V ai = - (V a + V b ) 



(20) 



Applying stability analysis to the resulting SPH equations, allows 
us to obtain a somewhat complicated dispersion relation. An 



appropriate form for the kernel can be derived in a similar fashion 
to that employed in the case of constant particle separation. 

10. Stability Analysis with Viscosity 

One might think that introducing viscosity to the equations 
of motion might eliminate the instability resulting from negative 
stress. In particular, since the first instability is of short wave- 
length (k % 7r/Ax) we would expect viscosity to have a strong 
damping effect upon it. Let us consider the previous equations 
with the addition of artificial viscosity, 



dv a 
dt 
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Once again, we consider a linear wave perturbing an infinite 
line of particles parallel to the a>axis and obtain the following 
dispersion relation, 

cu 2 + iacu - b = 0, (22) 

(23) 



where, 
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The definitions of a and b are chosen such that they are real and, 
for the standard kernel and dt = 1, are positive. Thus, we have, 



(46 - a 2 



(25) 



From Eq. (5) we And, 



. .Ub-a 2 ) 2 a 
V exp [ikx — iwt) = V exp [ ikx =p % 1 — —t 

. ' (26) 

We see that a is responsible for dampening the motion and 
reducing the speed of propagation. Let us now consider the 
crucial wavenumber, k = 7r/Ax. Substituting into Eqs. (23) and 
(24) we obtain, 
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Typically if dt < (ie: negative stress) then b < 0. Thus, we 
obtain, 

.UP + a 2 ) 1 ' .a , 

w = ±i±— '- i-, (28) 

2 2' V ; 

where /? = —b. The component of this which leads to the 
instability is, 

L o= l -({Af} + a 2 )- - a), (29) 

So we see that the instability is not removed by the introduction 
of viscosity. If a 2 /A(3 < 1 then 



8/3 J 



(30) 



If a/2 > a 2 /8j3, then we see that the effect of introducing a small 
amount of viscosity is to reduce the growth rate of the instability. 

11. Two Dimensional Stability Analysis 

Let's consider the stability properties of standard two di- 
mensional isothermal SPH, 

d Vab {Pab . Pcd\- w 

ai c d \Pab Pcd/ 

Pab = c 2 p ab , p ab = ^2^2 m cd W abcd ( 31 ) 

c d 

W abcd = W (x ab - x cd , y ab - y cd , h) 

We can analyse the stability of a rectangular lattice of particles 
by considering the propagation of plane waves on such a grid, 

x ab = aAx + X exp (ik x aAx + ik y bAy — iwt) , (32) 
y ab = bAy + Y exp (ik x aAx + ik y bAy — iwt) . 

Proceeding in a similar fashion to the one dimensional analysis we 
obtain two solutions for each pair of k x and k y . The full details 
are in Morris 1994, but let us consider the dispersion relation 
obtained for plane waves propagating along the a>axis. We find 
that there is a longitudinal wave solution, 

to 2 = 2 > > (1 - coskAxi) — — 

Po z -^ z -^ v ' dx 2 




which corresponds to the one dimensional result Eq. (6) with 
dt = 1 and the kernel replaced by a sum over j of the kernel. 
It turns out that this finite approximation to the y integration 
of the kernel introduces some small instabilities for large hj Ax. 
The system also supports a transverse wave with the dispersion 
relation, 

2 mc 2 , . d 2 W / , 
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These waves, of course, are not present in the exact solution to 
the wave equation for an isothermal gas, but they are present 
in our lattice approximating an isothermal gas. Their presence 
is the result of variations in the value of the gradient of the 
kernel as particles "jostle" around. It turns out (see Morris 1994) 
that for spline interpolated kernels w 2 is negative for about half 
of the choices of h for a given Ax. The growth rate of these 
instabilities is radically reduced as higher order spline interpolants 
are used. If the Gaussian is employed the growth rate of these 
unstable transverse modes is negligible for all practical choices 
of smoothing length. This would suggest that, once again, the 
stability properties of SPH are improved by using kernels whose 
Fourier transform falls off more rapidly. 

It is interesting to note that the standard implementation of 
viscosity has no effect on linear transverse waves. The stability 
properties of hexagonal lattices have also been considered in 
Morris 1994 with somewhat similar results. Three dimensional 
stability analysis also reveals unstable transverse modes which 
are reduced by the use of higher order spline interpolants. It can 
also be shown that the pressure differencing formulation does not 
exhibit these instabilities. 

12. Conclusion 

It is observed that the stability properties of SPH are im- 
proved by the use of kernels whose Fourier transforms fall off 
more rapidly. The instability first observed in exactly momentum 
conserving SPMHD in Phillips and Monaghan 1985 can be re- 
produced by a simple artificial equation of state (Eq. (4)). There 



are many applications which involve an equation of state having 
the effect of adjusting the background pressure acting between 
particles. Some examples are MHD, incompressible SPH, elastic- 
plastic now and problems where we model a fluid experiencing 
an external pressure. The pressure adjustment, when SPH is im- 
plemented in an exactly momentum conserving form, influences 
the propagation of waves within the medium but this influence 
is reduced by the use of higher order spline interpolated kernels. 
Once the stress acting between particles changes sign, however, 
a different kernel must be used. It is possible to construct a 
kernel which varies with the stress acting between particles so 
as to ensure stability and realistic wave propagation. However, 
for complicated equations of state in two or three dimensions, 
this may not be feasible. Formulations of SPH which calculate 
the stress gradients by taking differences between the stress at 
neighbouring particles are observed to simulate the propagation 
of waves very well. The simulation of strong shocks, however, 
suffers since momentum is no longer conserved exactly. In many 
circumstances it may be best to split the stress tensor into a com- 
ponent which is always positive and evaluate the remainder using 
a differencing formulation. For example, in the case of MHD the 
magnetic pressure could be added to the hydrodynamic pressure 
and employed in an exactly momentum conserving form. This 
would permit the simulation of strong shocks for which pressure 
forces are dominant. It may also be possible to adjust the stress 
at each of the particles by a constant, so as to ensure it is posi- 
tive everywhere, thus permitting exact momentum conservation. 
However, if the adjustment is extreme, a higher order spline inter- 
polated kernel my be required in order to avoid instabilities and 
strong dispersive effects. Boundary conditions may also become 
complicated. 

For all two dimensional and three dimensional applications 
using an exactly momentum conserving form, the use of kernels 
with compact support introduces instabilities in transverse modes 
on rectangular and hexagonal lattices of particles. The growth 
rates of these instabilities are observed to decrease dramatically as 
higher order spline approximations to the Gaussian are employed. 
These instabilities are not exhibited by differencing formulations. 

In general, kernels more closely approximating a Gaussian 
will give better results, but will cost more computationally as the 
number of contributing neighbours increases. However, this cost 
may be offset since such kernels may permit a decrease in h. It is 
important to have an understanding of the detail one wishes to 
resolve in a given problem. This understanding combined with 
sound knowledge of the stability and accuracy of the method 
employed allows us to have confidence in the results obtained. 
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